Large-scale Monte Carlo simulations of the three-dimensional XY spin glass 
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We study the XY spin glass by large-scale Monte Carlo simulations for sizes up to 24 , down 
to temperatures below the transition temperature found in earlier work. The data for the larger 
sizes show more marginal behavior than that for the smaller sizes indicating that the lower critical 
dimension is close to, and possibly equal to three. We find that the spins and chiralities behave in 
a similar manner. We also address the optimal ratio of "over-relaxation" to "Metropolis" sweeps in 
the simulation. 
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I. INTRODUCTION 

Following the convincing numerical work of Ballesteros 
et ali there has been little doubt that Ising spin glasses 
in three dimensions have a finite temperature transition. 
In this paper we shall study a related model for which the 
existence of a finite temperature transition is more con- 
troversial: the isotropic XY spin glass, which is composed 
of classical spins with two components. Early work on 
this model in three dimensions 2,3 indicated a zero tem- 
perature transition, or possibly a transition at a very 
low but non-zero temperature. However, following the 
pioneering work of Villain^, which emphasized the role 
of "chiralities" (Ising-like variables which describe the 
handedness of the non-collinear spin structures), Kawa- 
mura and Tanemura 5 proposed that the spin glass tran- 
sition only occurs at T§q = and that a chiral glass 
transition occurs at a finite temperature Tea- This sce- 
nario requires that spins and chiralities decouple at long 
length scales. Kawamura and collaborators have given 
numerical results in favor of this scenario^. 

However, the absence of a spin glass transition in the 
XY spin glass has been challenged by Maucourt and 
Grempel? and subsequently Akino and Kosterlitz^ who 
found evidence for a possible finite T$g from zero temper- 
ature domain wall calculations. Furthermore, by study- 
ing the dynamics of the XY spin glass in the phase rep- 
resentation, Granato^ found that the "current-voltage" 
characteristics exhibited scaling behavior which he inter- 
preted as a transition in the spins as well as the chirali- 
ties. 

In earlier workiS,, referred to as LY, Lee and one of the 
present authors studied spin and chiral correlations on 
an equal footing, using the method of analysis that was 
the most successful for the Ising spin glass^ i 11 ' 12 , namely 
finite-size scaling of the correlation length. Considering 
a modest range of sizes, N = L 3 with L < 12, LY found 
that the behavior of spins and chiralities was quite sim- 
ilar and they both had a finite temperature transition, 
apparently at the same temperature. 

LY studied both XY and Heisenberg models, finding 
similar conclusions for both. However, for the Heisenberg 



case, subsequent studies on much larger size: 



,13.14 



up to 



L = 32, have painted a more complex picture. The data 



at the lowest temperatures and largest sizes seems rather 
"marginal" , i.e. the system is close to the lower criti- 
cal dimension where the finite-temperature phase tran- 
sition is removed by fluctuations. The data for spins 
and chiralities are still quite similar, though not iden- 
tical, and do not seem to give compelling evidence for 
spin-chirality decoupling as proposed by Kawamura. In 
addition, Hukushima and Kawamura^ have also stud- 
ied somewhat larger sizes than LY (L < 20), but they 
argued that their data is consistent with spin-chirality 
decoupling. 



It is of interest to know whether the "crossover" to 
more marginal behavior found for larger sizes is special 
to the three-component case, or whether the same sit- 
uation occurs quite generally with vector spin glasses. 
In this paper, we therefore study the XY (2-component) 
spin glass for larger sizes (up to 24 3 ) than in LY (which 
went only up to 12 3 ). We find a situation that is quite 
similar to the Heisenberg case, namely marginal behavior 
for low-T and large sizes. The behavior of the spin glass 
and chiral glass correlation length is very similar, more 
similar than was the case for the Heisenberg spin glass, 
and docs not appear to provide evidence for spin-chirality 
decoupling, at least up to the sizes studied. 



Simulations on very large sizes for vector spin glasses 
have been possible because including "overrelaxation" 
moves, in addition to the more familiar Metropolis or 
heatbath moves, speeds up equilibration 1 ^. A second 
motivation of the present work is to investigate quanti- 
tatively the optimal ratio of overrelaxation to Metropolis 
sweeps for the XY spin glass. 



The layout of this paper is as follows. Section |TT] de- 
scribes the model, the parameters of the simulations, and 
the finite-size scaling approach. The results for the cor- 
relation length are presented in Sec. IIIIl In Sec. IIVI wc 
estimate the optimal ratio between the number of overre- 
laxation and Metropolis sweeps, and Sec. [V] summarizes 
our conclusions. 
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II. MODEL AND ANALYSIS 

We use the standard Edwards- Anderson XY spin glass 
model 



n 



Jr/Sj ■ s 

(i,3) 



.1 1 



(i) 



where the Si are 2-component classical vectors of unit 
length at the sites of a simple cubic lattice, and the Jij 
are nearest neighbor interactions with a Gaussian dis- 
tribution with zero mean and standard deviation unity. 
Periodic boundary conditions are applied on lattices with 
N = L 3 spins. 

The spin glass order parameter, ^"(k), at wave vector 
k, is defined to be 



(2) 



where fi and v are spin components, and "(1)" and "(2)" 
denote two identical copies of the system with the same 
interactions. From this we determine the wave vector 
dependent spin glass susceptibility xsgO*) by 



Xsg 



(k) = iv£[(|<r(k)f 



(3) 



where (• • •) denotes a thermal average and [• • -] av de- 
notes an average over disorder. The spin glass correlation 
length is then determined^^ from 



1 



Xsg(O) 
2sin(fc min /2) \ 

\ XSG (H-min ) 



- 1 



1/2 



(4) 



where k min = (2tt/£)(1, 0, 0). 

For the XY spin glass, chirality of a square is 6 

1 ' 

«f = — = ^ sgn( J /m ) sin(6» ; - 8 m ), 

(Z,m) 



(5) 



where £?/ is the angle characterizing the direction of spin 
Si , and the prime on the sum indicates that it is over the 
four bonds around the elementary plaquette perpendic- 
ular to the /x axis whose "bottom left" corner is site i. 
The chiral glass susceptibility is then given by 

X £ G (k)=A[(|^(k)| 2 )] av , (6) 

where the chiral overlap (k) is given by 

^(k) = ^E< (1) < (2)eik - R< - ( 7 ) 



We define the chiral correlation lengths ^ L by 



2sin(fc min /2) \ 



' 1 y /2 , (8) 



in which xcG(k = 0) is independent of /i. Note that 
^ L will, in general, be different for /t along k m j n (the x 
direction) and perpendicular to k, though this difference 
is very small for large sizes. The results presented will 
be an average over the three (two transverse and one 
longitudinal) correlation lengths. 

To equilibrate the system efficiently we perform three 
types of Monte Carlo move. 

Firstly we use "over-relaxation" sweeps^ in which we 
sweep sequentially through the lattice, and, at each site, 
compute the local field on the spin, = Y]- Sj . The 
new value for the spin on site i is taken to be its old value 
reflected about H, i.e. 



Si 



_ Sj ■ H; 

2 ., H, 



(9) 



Over-relaxation sweeps preserve energy and so are also 
known as microcanonical sweeps. 

Secondly, we include Metropolis sweeps since, unlike 
the over-relaxation sweeps, these do change the energy, 
and so are needed to bring the system to equilibrium. For 
the data presented in Sees . HT1 and IIIII we do one Metropo- 
lis sweep after every 10 over-relaxation sweeps. As for 
the over-relaxation case, we sweep sequentially through 
the lattice. To update a given spin, we choose a trial 
new direction randomly within a window ±A0/2 of the 
current direction, and accept this new direction with the 
usual Metropolis probability, min(l, exp(—f3AE)), where 
/3 = 1/T and AE is the energy difference between the 
trial state and the current state. We choose the window 
size A9 to vary with temperature in such a way that the 
acceptance ratio for Metropolis moves is in the range of 
30 to 50%. 

A Metropolis sweep requires more CPU time than a 
over-relaxation sweep, so we do mainly over-relaxation 
sweeps, including some Metropolis sweeps only to change 
the energy from time to time to ensure that the algorithm 
is ergodic. In fact, as discussed in Sec. HV| including a 
fraction of over-relaxation sweeps not only reduces the 
CPU time (for a given total number of sweeps) but also 
reduces the number of sweeps needed to equilibrate. 

Finally we do "parallel tempering" sweep s 18 ! 19 , which 
are necessary to prevent the system being trapped in a 
valley in configuration space at low temperatures. One 
takes Nt copies of the system with the same bonds but at 
a range of different temperatures. The minimum temper- 
ature, T m j n = T\ , is the low temperature where one wants 
to investigate the system (below Tsg m our case), and 
the maximum, T max = T^ T , is high enough that the the 
system equilibrates very fast (well above Tsg m our case). 
A parallel tempering sweep consists of swapping the tem- 
peratures of the spin configurations at a pair of neighbor- 
ing temperatures, Ti and TJ+i, for i = 1, 2, • • • , T/v T _i 
with a probability that satisfies the detailed balance con- 
dition. Further details on th e ap plication to vector spin 
glasses can be found in Ref. LLJ. For the simulations in 
SecslTTl and IIIII we do one parallel tempering sweep after 
each Metropolis sweep. 
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TABLE I: Parameters of the simulations described in Sees. [II] 
and Mil iVsamp is the number of samples, iVS^i is the num- 
ber of over-relaxation Monte Carlo sweeps for equilibration 
for each of the 2Nt replicas for a single sample, and iV^iL i s 
the number of over-relaxation sweeps for measurement. The 
number of Metropolis sweeps and the number of parallel tem- 
pering sweeps are both equal to 10% of the number of over- 
relaxation SW66pS. Tin in 3-11(1 Tmax 

are the lowest and highest 
temperatures simulated, and Nt is the number of tempera- 
tures used in the parallel tempering. 
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FIG. 1: (Color online) Equilibration plot, testing Eq. (flO|| . 
for L — 16 at T — 0.265. It is seen that the data for U and 
U(qi,q s ), given by Eq. (|11[) . come together when the total 



number of over-relaxation sweeps, iV s , 



A* 



OR , ArOfl 
cquil T J'meas 



see Table HI is equal to about 2 x 10 5 . These two quantities 
then stay at their common value indicating that equilibration 
has been achieved. It is seen that the energy comes close to 
its equilibrium value very quickly, whereas U(qi,q s ), which 
depends on the link overlap qi between two replicas, takes 
much longer. 



Table U gives the parameters of the simulations used to 
collect the data in Sees. [TTI and IIIII 

To test for equilibration 2 ^ we require that data satisfy 




N 0R 
sweep 



FIG. 2: (Color online) A plot of the spin glass and chiral glass 
correlation lengths, £l and £l,c, divided by L, as a function of 
the total number of sweeps for L = 16 at T — 0.265. It is seen 
the data flattens off at around 2 x 10 sweeps, the value where 
the two sets of data in Fig.[T]start to agree. This indicates that 
when the data in Fig. [1] agree within high precision, i.e. when 
Eq. (|10|l is satisfied, the correlation lengths have reached their 
equilibrium value. 



the relation^ 



where 



U = U(q h q a ) 



U{qi,q s ) = 2^ (<U -9s), 



(10) 



(11) 



which is valid for a Gaussian bond distribution. Here 



j is the average energy per spin, 
) 2 ] av is the "link overlap", q s = 
iV , N b = (z/2)N is the number 



u = -\E 

8 = (lM)E(y)[(Si'S ; 

of nearest neighbor bonds, and z (= 6 here) is the lattice 
coordination number. Equation (|10[) is easily derived by 
integrating by parts the expression for the average energy 
with respect to Jy, noting that the average [• • -] av is over 
a Gaussian function of the Jy's. 

The spins are initialized in random directions so the 
energy, the LHS of Eq. (fTOjl , is initially close to zero and 
decreases, presumably monotonically, to its equilibrium 
value as the length of the simulation increases. Hence the 
LHS of Eq. (fTU]) will be too large if the simulation is too 
short to equilibrate the system. On the other hand, the 
RHS of Eq. (HUJ), will be too small if the simulation is too 
short because qi starts off close to zero and then increases 
with MC time as the two replicas start to find the same 
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local minima. The quantity q s will be less dependent on 
Monte Carlo time than qi since it is a local variable for 
a single replica. (For the Ising case it is just a constant.) 
Hence if the simulation is too short the RHS of Eq. (fTU|) 
will be too low. In other words, the two sides of Eq. (fT0|) 
are expected to approach the common equilibrium value 
from opposite directions as the length of the simulation 
increases. Only if Eq. (|10[) is satisfied within small error 
bars do we accept the results of a simulation. 

Figure [1] shows a test to verify that Eq. (fT0|) is sat- 
isfied at long times. For the parameters used, L = 
16, T — 0.265, this occurs when the total number of (over- 
relaxation) sweeps (Ng^ cp = N^ n + N^) is about 
2 x 10 5 . Figure [5] shows that the spin and chiral corre- 
lation lengths appear to become independent of N SVfeep , 
and hence are presumably equilibrated, when N svreep is 
larger than this same value. Hence, it appears that when 
Eq. (flU|) is satisfied to high precision, the data for the 
correlation lengths is equilibrated. 

With the number of sweeps shown in Table QJ Eq. (flQ|) 
was satisfied for all sizes and temperatures. The error 
bars are made sufficiently small by averaging over a large 
number of samples. 

Since £l /L is dimensionless it has the finite size scaling 
formal 
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x(l x i v {t-t sg ) 



(12) 
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FIG. 3: (Color online) Data for £l/L, the spin glass corre- 
lation length divided by system size, as a function of T for 
different system sizes. 



where v is the correlation length exponent. Note that 
there is no power of L multiplying the scaling function 
X. By contrast, for the spin glass susceptibility, \sg = 
XSG(k = 0), which has dimensions, the finite-size scaling 
form is 



XSG = l 2 -" sb K (l}/ v {T - T SG ) 



(13) 



where rjsc is a critical exponent. There is an expression 
analogous to Eq. (|12p for the chiral correlation length, 
and to Eq. lfT5]) for the chiral glass susceptibility xcg = 
Xccfe = 0)- For the later case, there is no reason to 
expect that the exponents tjsg and ?7cg are equal. 

From Eq. (fl~2)) it follows that the data for £l/L for 
different sizes come together at T = Tqq. In addition, 
they are also expected to splay out again on the low-T 
side 1 if there is spin glass order below T$g- I n a marginal 
situation with a line of critical points, as in the Kosterlitz- 
Thouless-Berezinskii theory of the transition in the two- 
dimensional XY ferromagnet, the data for different sizes 
would come together at Tsg and then stick together at 
lower T, see for example Fig. 3 of Ref. [l|. 
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III. RESULTS 

We studied sizes from L = 4 to L = 24, as shown in 
Table HI The CPU time involved to get this data is about 
8 Mac G5 CPU years. 



FIG. 4: (Color online) Data for the chiral correlation length 
(averaged over longitudinal and transverse directions) divided 
by system size, as a function of T for different system sizes. 
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FIG. 5: (Color online) The same data as in Figs. [3] and [4] but 
including only the largest sizes and in a somewhat expanded 
scale. 



The data for the spin glass correlation length (divided 
by L) is shown in Fig. [31 and the corresponding data for 
the chiral glass correlation length is shown in Fig. [U In 
both cases the data for smaller sizes intersect and splay 
out at lower temperature. However, for the larger sizes 
the splaying out is small, indicating close to "marginal" 
behavior, i.e. the "lower critical dimension" is close to 3. 

The data for the spins and chiralities in Figs. [3] and 
Fig. 01 are very similar, so we do not see evidence for 
spin-chirality decoupling. To make clearer the similarity 
between the two sets of data we plot them both in Fig.[5l 
including just the three largest sizes. The temperature 
where the data merges decreases slightly with increasing 
size. We have estimated the temperatures where the data 
intersect/merge for different pairs of sizes and present 
the results in Table [TTJ The temperatures are seen to 
decrease with increasing size. If one neglects the smallest 
pair of sizes (L = 4/6) the shift is somewhat bigger for 
the spins than for the chiralities but, from the data, it 
is not possible to reliably estimate whether or not the 
intersection temperature will tend to zero for L — > oo for 
either set of data. 



TABLE II: Estimated crossing temperatures for the spin and 
chiral glass correlation lengths. The results are given to the 
nearest 0.005, but the uncertainties are greater than this be- 
cause of the error bars in the data itself. 



sizes Tci-ossing (spins) Tossing (chiralities) 

4/6 0.355 0.375 

6/8 0.33 0.32 

8/12 0.33 0.335 

12/16 0.31 0.32 

16/24 0.285 O30 



In Fig. [6] we present data for the ratio of the chiral- glass 
to spin-glass correlation lengths. For the largest sizes the 
data intersects for T about 0.33 and then (slightly) splays 
out in the low-T side. If there is a single transition in- 
volving both spins and chiralities, then the data would 
become independent of size at the transition (since both 
£l and £ Ci l are proportional to L there, see Eq. ([T2]) ). If 
the stiffness exponents for spins and chiralities are equal 
(we are not aware of any argument for this even if there is 
a single transition) then the data would become indepen- 
dent of L for large L at low-T. If the stiffness exponent 
for chiralities is larger than that for the spins, then the 
ratio would diverge in this limit. From the data it is 
not possible to say for sure if the data diverges or not at 
low-T, but the size dependence at the larger sizes is very 
weak. 

In the spin-chirality decoupling scenario, the ratio 
would diverge even at the transition, and there would 
not be a common intersection. We feel that the data of 
Fig. [6] reinforces our view that if spin-chirality decoupling 
occurs one would need even larger sizes than L — 24 to 
see it. 

We also present data for the spin-glass and chiral-glass 
susceptibilities in Figs. [7] and [5] respectively. Dividing 
by L 2 ~'', where r/ is a critical exponent, the data should 
intersect at the critical temperature, see Eq. (IT31) . where 
■qsa is not necessarily equal to t\cg- ln order to get 
intersections for T ~ 0.30, where the correlation data 
merge/intersect for the largest sizes, we took -qsa = —0.2 
and tjcg = 0.1 in the plots. 

Given the large corrections to scaling clearly visible in 
the data for the correlation lengths, it does not appear 
possible to get reliable estimate of the critical exponents, 
rjsa and Vcc, or of the correlation length exponent v. 



IV. OPTIMIZING THE FRACTION OF 
OVERRELAXATION SWEEPS 

As already noted, adding overrelaxation steps has been 
observe d 13 ' 14 i 16 to speed up equilibration. Here we look 
systematically at how the ratio of the number of over- 
relaxation (OR) sweeps to Metropolis (MET) sweeps al- 
ters the total number of sweeps needed to equilibrate. 
In Fig. [HI we plot both sides of Eq. (fTU)) . which are 
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FIG. 6: (Color online) Data for the ratio of the chiral glass 
to the spin glass correlation lengths for sizes from 8 to 24. 




FIG. 7: (Color online) Data for the spin glass susceptibility 
Xsg = XSG(k = 0) divided by L 2 ~ VSG where we took tjsg = 
—0.2 in order to get the data to intersect (see Eq. (|13p ) for 
T around 0.30 since this is roughly where the data for £l/L 
and £ c l intersect/merge for the largest sizes, see Figs. [3] and 

H 



FIG. 8: (Color online) Similar to Fig.Obut for the chiral glass 
susceptibility xcg = XCG(k = 0). Here we took tjcg = 0.1. 



equal in equilibrium, for different ratios of the num- 
ber of OR sweeps to MET sweeps. The data is for 
L = 16, T = 0.265. It is seen that equilibration is con- 
siderably speeded up by including OR sweeps. It seems 
that doing 10 OR per MET (which was used in the re- 
sults in the earlier sections) is somewhat better than 1 
OR or 40 OR. Reference [Tj argues that of order L OR 
sweeps should be done for each MET sweep "to let the 
microcanonical wave run over the system" . Our data is 
consistent with this, though it seems that time to equili- 
brate is not very sensitive to the precise ratio of OR to 
MET sweeps. 

We should emphasize that including OR sweeps not 
only reduces the number of sweeps to equilibrate, as seen 
in Fig. [31 but also reduces the CPU time by an even 
bigger factor, because each OR sweep runs several times 
faster on the computer than an MET sweep. 



V. CONCLUSIONS 

We have studied the XY spin glass in three dimen- 
sions by Monte Carlo simulations using larger sizes than 
before. We find that the lower critical dimension is close 
to three. We also find that the behavior of the spin glass 
and chiral glass correlations lengths is strikingly simi- 
lar, see Fig. [5J and, in our view, does not support the 
spin chirality decoupling scenario, at least for sizes up to 
L = 24. 

In earlier work, Maucourt and Grempcl 7 have studied 
the 3d XY spin glass using the domain- wall renormaliza- 
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FIG. 9: (Color online) Results for L = 16, T = 0.265. The 
data connected by solid lines is U(qi,q s ) in Eq. (|ll|l for dif- 
ferent number of over-relaxation (OR) sweeps per Metropolis 
(MET) sweep as indicated. The horizontal axis is the total 
number of OR plus MET sweeps. The data connected by the 
dashed line is the energy U, which should equal U(qi,q 3 ) in 
equilibrium according to Eq. (|10[> . Since the energy equili- 
brates relatively fast, its value does not depend significantly 
on the ratio of OR to MET sweeps for the range of sweeps 
presented. The number of parallel tempering sweeps is the 
same for all sets of data except for "40 OR" where it is 1/4 
as many. 

tion group (DWRG), for sizes up to L = 8. They argue 
that there is a positive stiffness for the chiralities, and 
hence a finite temperature transition, while for spin glass 
ordering the system is close to its lower critical dimen- 
sion. The conclusion for chiralities is different from ours 
but we note that our sizes are much larger (L < 24) and 
that we only see marginal behavior in the chiralities for 
L > 12. Furthermore, our approach gives directly the 
correlation lengths, whereas for the DWRG ground state 
energies with different boundary conditions are computed 
from which a stiffness is inferred. 

Kawamura and Li 6 used Monte Carlo simulations with 



sizes up to L = 16 to compute the overlap function of 
the spins and chiralities. In particular, they compute 
the "Binder ratio" which, like the ratio of the correla- 
tion length to system size studied here, is dimensionless. 
The spin glass Binder ratio is found to monotonically de- 
crease with increasing L at each temperature. However, 
we feel that use of the Binder ratio can be tricky near 
the lower critical dimension especially when the num- 
ber of components of the order parameter is high. Since 
the spin glass order parameter is quadratic in the spins 
and the spins have two components, the order param- 
eter has four-components here. The Binder ratio looks 
at the change in shape of the distribution of the (square 
root of the) order parameter squared summed over all 
components, when going below the transition. Because 
of the central limit theorem, there would be no change 
in shape for an infinite number of components. If the 
number is large the change in shape is small and can eas- 
ily be masked by corrections to scaling, especially if the 
system is close to the lower critical dimension where cor- 
rections only fall off very slowly with system size. The 
use of the Binder ratio for vector spin glasses has also 
been criticized by Shirakura and Matsubara^i (they con- 
sidered explicitly the Heisenberg case). For the chiral 
glass Binder ratio, Kawamura and Li estimate a transi- 
tion temperature from a dip in the data. However, even 
if the transition is of an unconventional kind (as they 
claim in order to explain the dip) it seems to us that the 
Binder ratio should still increase with increasing L at low 
temperature if there is chiral glass order. However, this 
is not observed. 

We therefore argue that our results, which compute di- 
rectly the relevant correlation lengths, indicate that spin- 
chirality decoupling does not seem to occur, at least for 
sizes up to L = 24. 

Finally, we find that equilibration is considerably 
speeded up by performing several (perhaps of order L) 
over-relaxation sweeps per Metropolis sweep, see Fig. GO 
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